package com.ColonyCount.Image;

import java.security.InvalidParameterException;
import java.util.Arrays;


import android.util.Log;

public class Segmentation {

	/** the following constants are used to set bits corresponding to pixel types */
	final static byte MAXIMUM = (byte)1;			// marks local maxima (irrespective of noise tolerance)
	final static byte LISTED = (byte)2;             // marks points currently in the list
    final static byte PROCESSED = (byte)4;          // marks points processed previously
    final static byte MAX_AREA = (byte)8;           // marks areas near a maximum, within the tolerance
    final static byte EQUAL = (byte)16;             // marks contigous maximum points of equal level
    final static byte MAX_POINT = (byte)32;         // marks a single point standing for a maximum
    final static byte ELIMINATED = (byte)64;        // marks maxima that have been eliminated before watershed
    
    /** what type of output to create was chosen in the dialog (see constants below)*/
	
	public final static int SINGLE_POINTS=0;		/** Output type single points */
	public final static int IN_TOLERANCE=1;			/** Output type all points around the maximum within the tolerance */
	public final static int SEGMENTED=2;			/** Output type watershed-segmented image */
	public final static int POINT_SELECTION=3;		/** Do not create image, only mark points */
	public final static int COUNT=4;				/** Do not create an image, just count maxima and add count to Results table*/
	int outputType = SEGMENTED;
    
	final static int IS_LINE=1;                     // a point on a line (as a return type of isLineOrDot)
    final static int IS_DOT=2;                      // an isolated point (as a return type of isLineOrDot)	
	
	int[] dirOffset;
	int[] dirXoffset;
	int[] dirYoffset;
	
	boolean excludeEdgesNow = false;
	boolean excludeOnEdges = true;
	
    float globalMin = Float.MAX_VALUE;
    float globalMax = -Float.MAX_VALUE;
	MaxPoint[] maxPoints;
		
	   /** returns whether the neighbor in a given direction is within the image
     * NOTE: it is assumed that the pixel x,y itself is within the image!
     * Uses class variables width, height: dimensions of the image
     * @param x         x-coordinate of the pixel that has a neighbor in the given direction
     * @param y         y-coordinate of the pixel that has a neighbor in the given direction
     * @param direction the direction from the pixel towards the neighbor (see makeDirectionOffsets)
     * @return          true if the neighbor is within the image (provided that x, y is within)
     */
	boolean isWithin(Image16 inputMap, int x, int y, int direction) 
	{
        int width = inputMap.getWidth();
        int height = inputMap.getHeight();
        int xmax = width - 1;
        int ymax = height -1;
        switch(direction) {
            case 0:
                return (y>0);
            case 1:
                return (x<xmax && y>0);
            case 2:
                return (x<xmax);
            case 3:
                return (x<xmax && y<ymax);
            case 4:
                return (y<ymax);
            case 5:
                return (x>0 && y<ymax);
            case 6:
                return (x>0);
            case 7:
                return (x>0 && y>0);
        }
        return false;   
    } // isWithin	
	
	int trueEdmHeight(int x, int y, Image16 inputMap) {
        int xmax = inputMap.getWidth()-1;
        int ymax = inputMap.getHeight()-1;
        int v =  inputMap.getPixelValue(x, y);
        if (x==0 || y==0 || x==xmax || y==ymax || v==0) {
            return v;                               //don't recalculate for edge pixels or background
        } else {
            int one = EDM.ONE;
            int sqrt2 = EDM.SQRT2;
            int trueH = v+sqrt2/2;                  //true height can never by higher than this
            boolean ridgeOrMax = false;
            for (int d=0; d<4; d++) {               //for all directions halfway around:
                int d2 = (d+4)%8;                   //get the opposite direction and neighbors
                int v1 = inputMap.getPixelValue(x+dirXoffset[d], y+dirYoffset[d]);
                int v2 = inputMap.getPixelValue(x+dirXoffset[d2], y+dirYoffset[d2]);
                int h;
                if (v>=v1 && v>=v2) {
                    ridgeOrMax = true;
                    h = (v1 + v2)/2;
                } else {
                    h = Math.min(v1, v2);
                }
                h += (d%2==0) ? one : sqrt2;        //in diagonal directions, distance is sqrt2
                if (trueH > h) trueH = h;
            }
            if (!ridgeOrMax) trueH = v;
            return trueH;
        }
    }
	
	
	class MaxPoint implements Comparable {
        float value;
        short x;
        short y;

        /** a constructor filling in the data */
        MaxPoint(short x, short y, float value) 
        {
            this.x = x;
            this.y = y;
            this.value = value;
        }

        /** a comparator required for sorting (interface Comparable) */
		 public int compareTo(Object o) 
		 {
		    //return Float.compare(value,((MaxPoint)o).value); //not possible since this requires Java 1.4
		    float diff = value-((MaxPoint)o).value;
		    if (diff > 0f) return 1;
		    else if (diff == 0f) return 0;
		    else return -1;
		}
    } 
	
	private void makeDirectionOffsets(int width, int height) 
	{
	    dirOffset = new int[] { -width, -width+1, +1, +width+1, +width, +width-1,   -1, -width-1 };
	    dirXoffset = new int[] {    0,      1,     1,     1,        0,     -1,      -1,    -1    };
	    dirYoffset = new int[] {   -1,     -1,     0,     1,        1,      1,       0,    -1,   };
    }
	  
	 MaxPoint[] getSortedMaxPoints(Image16 inputMap, Image8 maxMap, boolean excludeEdgesNow, boolean isEDM, float globalMin, double threshold) 
	 {
         int width = inputMap.getWidth();
         int height = inputMap.getHeight();

         byte[] types =  maxMap.pixels;
         int nMax = 0;  //counts local maxima

         boolean checkThreshold = true;
         
         for (int y=0, i=0; y<height; y++) 
         {         // find local maxima now
             for (int x=0; x<width; x++, i++) 
             {      // for better performance with rois, restrict search to roi
                 float v = inputMap.getPixelValue(x, y);
                 float vTrue = isEDM ? trueEdmHeight(x, y, inputMap) : v;  // for 16-bit EDMs, use interpolated ridge height
                 if (v==globalMin) continue;
                 if (excludeEdgesNow && (x==0 || x==width-1 || y==0 || y==height-1)) continue;
                 if (checkThreshold && v<threshold) continue;
                 boolean isMax = true;
                 /* check wheter we have a local maximum.
                  Note: For a 16-bit EDM, we need all maxima: those of the EDM-corrected values
                  (needed by findMaxima) and those of the raw values (needed by cleanupMaxima) */
                 for (int d=0; d<8; d++) 
                 {                         // compare with the 8 neighbor pixels
                     if (isWithin(inputMap, x, y, d)) {
                         float vNeighbor = inputMap.getPixelValue(x+dirXoffset[d], y+dirYoffset[d]);
                         float vNeighborTrue = isEDM ? trueEdmHeight(x+dirXoffset[d], y+dirYoffset[d], inputMap) : vNeighbor;
                         if (vNeighbor > v && vNeighborTrue > vTrue) {
                             isMax = false;
                             break;
                         }
                     }
                 }
                 if (isMax) {
                     types[i] = MAXIMUM;
                     nMax++;
                 }
             } // for x
         } // for y
         
         MaxPoint [] maxPoints = new MaxPoint[nMax];
         int iMax = 0;
         for (int y=0, i=0; y < height; y++)           //enter all maxima into an array
             for (int x=0; x < width; x++, i++)
                 if (types[i]==MAXIMUM) {
                     maxPoints[iMax] = new MaxPoint((short)x, (short)y, isEDM?trueEdmHeight(x,y,inputMap):inputMap.getPixelValue(x,y));
                     //if(ip.getPixelValue(x,y)>=250) IJ.write("x,y,value="+maxPoints[iMax].x+","+maxPoints[iMax].y+","+maxPoints[iMax].value);
                     iMax++;
                 }
         Arrays.sort(maxPoints);                                 //sort the maxima by height
         return maxPoints;
     } 
	
	 /** Check all maxima in list maxPoints, mark type of the points in typeP
	    * @param ip             the image to be analyzed
	    * @param typeP          8-bit image, here the point types are marked by type: MAX_POINT, etc.
	    * @param maxPoints      input: a list of all local maxima, sorted by height
	    * @param excludeEdgesNow whether to avoid edge maxima
	    * @param isEDM          whether ip is a 16-bit Euclidian distance map
	    * @param globalMin      minimum pixel value in ip
	    * @param tolerance      minimum pixel value difference for two separate maxima
	    */   
	 void analyzeAndMarkMaxima(Image16 inputMap, Image8 maxMap, MaxPoint[] maxPoints, boolean excludeEdgesNow, boolean isEDM, float globalMin, double tolerance) 
	 {
		 long time = System.currentTimeMillis();
		 
        int width = inputMap.getWidth();
        int height = inputMap.getHeight();
        byte[] types =  maxMap.pixels;
        int nMax = maxPoints.length;
        short[] xList = new short[width*height];    //here we enter points starting from a maximum
        short[] yList = new short[width*height];
//        Vector xyVector = null;
//        Roi roi = null;
        boolean displayOrCount = (outputType==POINT_SELECTION||outputType==COUNT);
        
        for (int iMax=nMax-1; iMax>=0; iMax--) 
        {    //process all maxima now, starting from the highest
            if (iMax%100 == 0 && Thread.currentThread().isInterrupted()) 
            	return;
            float v = maxPoints[iMax].value;
            if (v==globalMin) 
            	break;
            int offset = maxPoints[iMax].x + width*maxPoints[iMax].y;
            if ((types[offset]&PROCESSED)!=0)       //this maximum has been reached from another one, skip it
                continue;
            xList[0] = maxPoints[iMax].x;           //we create a list of connected points and start the list
            yList[0] = maxPoints[iMax].y;           //  at the current maximum
            types[offset] |= (EQUAL|LISTED);        //mark first point as equal height (to itself) and listed
            int listLen = 1;                        //number of elements in the list
            int listI = 0;                          //index of current element in the list
            boolean isEdgeMaximum = (xList[0]==0 || xList[0]==width-1 || yList[0]==0 || yList[0]==height-1);
            boolean maxPossible = true;             //it may be a true maxmum
            double xEqual = xList[0];               //for creating a single point: determine average over the
            double yEqual = yList[0];               //  coordinates of contiguous equal-height points
            int nEqual = 1;                         //counts xEqual/yEqual points that we use for averaging
            do {
                offset = xList[listI] + width*yList[listI];
                for (int d=0; d<8; d++) {           //analyze all neighbors (in 8 directions) at the same level
                    int offset2 = offset+dirOffset[d];
                    if (isWithin(inputMap, xList[listI], yList[listI], d) && (types[offset2]&LISTED)==0) {
                        if ((types[offset2]&PROCESSED)!=0) {
                            maxPossible = false;    //we have reached a point processed previously, thus it is no maximum now
                            //if(xList[0]>510&&xList[0]<77)IJ.write("stop at processed neighbor from x,y="+xList[listI]+","+yList[listI]+", dir="+d);
                            break;
                        }
                        int x2 = xList[listI]+dirXoffset[d];
                        int y2 = yList[listI]+dirYoffset[d];
                        float v2 = inputMap.getPixelValue(x2, y2);
                        if (isEDM && (v2 <=v-(float)tolerance)) v2 = trueEdmHeight(x2, y2, inputMap); //correcting for EDM peaks may move the point up
                        if (v2 > v) {
                            maxPossible = false;    //we have reached a higher point, thus it is no maximum
                            //if(xList[0]>510&&xList[0]<77)IJ.write("stop at higher neighbor from x,y="+xList[listI]+","+yList[listI]+", dir,value,value2,v2-v="+d+","+v+","+v2+","+(v2-v));
                            break;
                        } else if (v2 >= v-(float)tolerance) {
                            xList[listLen] = (short)(x2);
                            yList[listLen] = (short)(y2);
                            listLen++;              //we have found a new point within the tolerance
                            types[offset2] |= LISTED;
                            if (x2==0 || x2==width-1 || y2==0 || y2==height-1) {
                                                isEdgeMaximum = true;
                                                if (excludeEdgesNow) {
                                    maxPossible = false;
                                          break;          //we have an edge maximum;
                                }
                            }
                            if (v2==v) {            //prepare finding center of equal points (in case single point needed)
                                types[offset2] |= EQUAL;
                                xEqual += x2;
                                yEqual += y2;
                                nEqual ++;
                            }
                        }
                    } // if isWithin & not LISTED
                } // for directions d
                listI++;
            } while (listI < listLen);
            byte resetMask = (byte)~(maxPossible?LISTED:(LISTED|EQUAL));
            xEqual /= nEqual;
            yEqual /= nEqual;
            double minDist2 = 1e20;
            int nearestI = 0;
            for (listI=0; listI<listLen; listI++) {
                offset = xList[listI] + width*yList[listI];
                types[offset] &= resetMask;     //reset attributes no longer needed
                types[offset] |= PROCESSED;     //mark as processed
                if (maxPossible) {
                    types[offset] |= MAX_AREA;
                    if ((types[offset]&EQUAL)!=0) {
                        double dist2 = (xEqual-xList[listI])*(xEqual-xList[listI]) + (yEqual-yList[listI])*(yEqual-yList[listI]);
                        if (dist2 < minDist2) {
                            minDist2 = dist2;   //this could be the best "single maximum" point
                            nearestI = listI;
                        }
                    }
                }
            } // for listI
            if (maxPossible) 
            {
                types[xList[nearestI] + width*yList[nearestI]] |= MAX_POINT;
                
//                if (displayOrCount && !(this.excludeOnEdges && isEdgeMaximum)) 
//                {
//                    if (xyVector==null) {
//                        xyVector = new Vector();
//                        roi = imp.getRoi();
//                    }
//                    short mpx = xList[nearestI];
//                    short mpy = yList[nearestI];
//                    if (roi==null || roi.contains(mpx, mpy))
//                        xyVector.addElement(new MaxPoint(mpx, mpy, 0f));
//                }
            }
        } 
        
//        if (displayOrCount && xyVector!=null) {
//            int npoints = xyVector.size();
//            if (outputType == POINT_SELECTION) {
//                int[] xpoints = new int[npoints];
//                int[] ypoints = new int[npoints];
//                MaxPoint mp;
//                for (int i=0; i<npoints; i++) {
//                    mp = (MaxPoint)xyVector.elementAt(i);
//                    xpoints[i] = mp.x;
//                    ypoints[i] = mp.y;
//                }
//                imp.setRoi(new PointRoi(xpoints, ypoints, npoints));
//            }
//            if (outputType==COUNT) {
//                ResultsTable rt = ResultsTable.getResultsTable();
//                rt.incrementCounter();
//                rt.setValue("Count", rt.getCounter()-1, npoints);
//                rt.show("Results");
//            }
//        }
        
        Log.v("", "Segmentation.analyzeAndMarkMaxima " + (System.currentTimeMillis() - time)+ "ms");
    } //void analyzeAndMarkMaxima	 
	 
	 static public void loadARGBFromTypeMap(Image8 typeMap, Image32 image)
	 {
		if(typeMap.getWidth() != image.getWidth() || typeMap.getHeight() != image.getHeight())
		{
			throw new InvalidParameterException("source and target not same size");
		}
		
		byte[] typePixels = typeMap.pixels;
		int[] imagePixels = image.pixels;
		int count = typePixels.length;

		for (int i = 0; i < count; i++) {
			int iVal = typePixels[i];
			
			if((iVal & ELIMINATED) > 0)
			{
				imagePixels[i] = 0x66FF0000;
				continue;
			}
			
			if((iVal & MAX_POINT) > 0)
			{
				imagePixels[i] = 0xFFFFFFFF;
				continue;
			}
			
			if((iVal & MAX_AREA) > 0)
			{
				imagePixels[i] = 0xFF999999;
				continue;
			}			
			
//			switch(iVal)
//			{
//				case 0: imagePixels[i] = 0x66000000;	
//				case MAXIMUM: imagePixels[i] = 0xFFFFFFFF; break;
//				case MAX_AREA: imagePixels[i] = 0xFF999999; break;
//				case EQUAL: imagePixels[i] = 0x660000FF; break;
//				case MAX_POINT: imagePixels[i] = 0x66FFFFFF; break;
//				default: imagePixels[i] = 0xFFFF0000; break;
//			}
//			
//			final static byte MAXIMUM = (byte)1;			// marks local maxima (irrespective of noise tolerance)
//			final static byte LISTED = (byte)2;             // marks points currently in the list
//		    final static byte PROCESSED = (byte)4;          // marks points processed previously
//		    final static byte MAX_AREA = (byte)8;           // marks areas near a maximum, within the tolerance
//		    final static byte EQUAL = (byte)16;             // marks contigous maximum points of equal level
//		    final static byte MAX_POINT = (byte)32;         // marks a single point standing for a maximum
//		    final static byte ELIMINATED = (byte)64;        // marks maxima that have been eliminated before watershed
					
			imagePixels[i] = 0xFF000000;
		}
	 }
	 
	    /** returns whether the neighbor in a given direction is within the image
	     * NOTE: it is assumed that the pixel x,y itself is within the image!
	     * Uses class variables width, height: dimensions of the image
	     * @param x         x-coordinate of the pixel that has a neighbor in the given direction
	     * @param y         y-coordinate of the pixel that has a neighbor in the given direction
	     * @param direction the direction from the pixel towards the neighbor (see makeDirectionOffsets)
	     * @return          true if the neighbor is within the image (provided that x, y is within)
	     */
	 boolean isWithin(Image8 ip, int x, int y, int direction) {
	        int width = ip.getWidth();
	        int height = ip.getHeight();
	        int xmax = width - 1;
	        int ymax = height -1;
	        switch(direction) {
	            case 0:
	                return (y>0);
	            case 1:
	                return (x<xmax && y>0);
	            case 2:
	                return (x<xmax);
	            case 3:
	                return (x<xmax && y<ymax);
	            case 4:
	                return (y<ymax);
	            case 5:
	                return (x>0 && y<ymax);
	            case 6:
	                return (x>0);
	            case 7:
	                return (x>0 && y>0);
	        }
	        return false;   //to make the compiler happy :-)
	    } // isWithin	 
	 
	 /** eliminate unmarked maxima for use by watershed. Starting from each previous maximum,
	     * explore the surrounding down to successively lower levels until a marked maximum is
	     * touched (or the plateau of a previously eliminated maximum leads to a marked maximum).
	     * Then set all the points above this value to this value
	     * @param outIp     the image containing the pixel values
	     * @param typeP     the types of the pixels are marked here
	     * @param maxPoints array containing the coordinates of all maxima that might be relevant
	     */    
	 void cleanupMaxima(Image8 outIp, Image8 typeP, MaxPoint[] maxPoints) {
	        byte[] pixels = (byte[])outIp.getPixels();
	        byte[] types = (byte[])typeP.getPixels();
	        int width = outIp.getWidth();
	        int height = outIp.getHeight();
	        int nMax = maxPoints.length;
	        short[] xList = new short[width*height];
	        short[] yList = new short[width*height];
	        for (int iMax = nMax-1; iMax>=0; iMax--) {
	            int offset = maxPoints[iMax].x + width*maxPoints[iMax].y;
	            //if (maxPoints[iMax].x==122) IJ.write("max#"+iMax+" at x,y="+maxPoints[iMax].x+","+maxPoints[iMax].y+": type="+types[offset]);
	            if ((types[offset]&(MAX_AREA|ELIMINATED))!=0) continue;
	            int level = pixels[offset]&255;
	            int loLevel = level+1;
	            xList[0] = maxPoints[iMax].x;           //we start the list at the current maximum
	            yList[0] = maxPoints[iMax].y;
	            //if (xList[0]==122) IJ.write("max#"+iMax+" at x,y="+xList[0]+","+yList[0]+"; level="+level);
	            types[offset] |= LISTED;                //mark first point as listed
	            int listLen = 1;                        //number of elements in the list
	            int lastLen = 1;
	            int listI = 0;                          //index of current element in the list
	            boolean saddleFound = false;
	            while (!saddleFound && loLevel >0) {
	                loLevel--;
	                lastLen = listLen;                  //remember end of list for previous level
	                listI = 0;                          //in each level, start analyzing the neighbors of all pixels
	                do {                                //for all pixels listed so far
	                    offset = xList[listI] + width*yList[listI];
	                    for (int d=0; d<8; d++) {       //analyze all neighbors (in 8 directions) at the same level
	                        int offset2 = offset+dirOffset[d];
	                        if (isWithin(typeP, xList[listI], yList[listI], d) && (types[offset2]&LISTED)==0) {
	                            if ((types[offset2]&MAX_AREA)!=0 || (((types[offset2]&ELIMINATED)!=0) && (pixels[offset2]&255)>=loLevel)) {
	                                saddleFound = true; //we have reached a point touching a "true" maximum...
	                                //if (xList[0]==122) IJ.write("saddle found at level="+loLevel+"; x,y="+xList[listI]+","+yList[listI]+", dir="+d);
	                                break;              //...or a level not lower, but touching a "true" maximum
	                            } else if ((pixels[offset2]&255)>=loLevel && (types[offset2]&ELIMINATED)==0) {
	                                xList[listLen] = (short)(xList[listI]+dirXoffset[d]);
	                                yList[listLen] = (short)(yList[listI]+dirYoffset[d]);
	                                listLen++;          //we have found a new point to be processed
	                                types[offset2] |= LISTED;
	                            }
	                        } // if isWithin & not LISTED
	                    } // for directions d
	                    if (saddleFound) break;         //no reason to search any further
	                    listI++;
	                } while (listI < listLen);
	            } // while !levelFound && loLevel>=0
	            for (listI=0; listI<listLen; listI++)   //reset attribute since we may come to this place again
	                types[xList[listI] + width*yList[listI]] &= ~LISTED;
	            for (listI=0; listI<lastLen; listI++) { //for all points higher than the level of the saddle point
	                offset = xList[listI] + width*yList[listI];
	                pixels[offset] = (byte)loLevel;     //set pixel value to the level of the saddle point
	                types[offset] |= ELIMINATED;        //mark as processed: there can't be a local maximum in this area
	            }
	        } // for all maxima iMax
	    } // void cleanupMaxima	 
	
	   /** Creates the lookup table used by the watershed function for dilating the particles.
	     * The algorithm allows dilation in both straight and diagonal directions.
	     * There is an entry in the table for each possible 3x3 neighborhood:
	     *          x-1          x          x+1
	     *  y-1    128            1          2
	     *  y       64     pxl_unset_yet     4
	     *  y+1     32           16          8
	     * (to find throws entry, sum up the numbers of the neighboring pixels set; e.g.
	     * entry 6=2+4 if only the pixels (x,y-1) and (x+1, y-1) are set.
	     * A pixel is added on the 1st pass if bit 0 (2^0 = 1) is set,
	     * on the 2nd pass if bit 1 (2^1 = 2) is set, etc.
	     * pass gives the direction of rotation, with 0 = to top left (x--,y--), 1 to top,
	     * and clockwise up to 7 = to the left (x--).
	     * E.g. 4 = add on 3rd pass, 3 = add on either 1st or 2nd pass.
	     */
	 private int[] makeFateTable() {
	        int[] table = new int[256];
	        boolean[] isSet = new boolean[8];
	        for (int item=0; item<256; item++) {        //dissect into pixels
	            for (int i=0, mask=1; i<8; i++) {
	                isSet[i] = (item&mask)==mask;
	                mask*=2;
	            }
	            for (int i=0, mask=1; i<8; i++) {       //we dilate in the direction opposite to the direction of the existing neighbors
	                if (isSet[(i+4)%8]) table[item] |= mask;
	                mask*=2;
	            }
	            for (int i=0; i<8; i+=2)                //if side pixels are set, for counting transitions it is as good as if the adjacent edges were also set
	                if (isSet[i]) {
	                    isSet[(i+1)%8] = true;
	                    isSet[(i+7)%8] = true;
	                }
	            int transitions=0;
	            for (int i=0, mask=1; i<8; i++) {
	                if (isSet[i] != isSet[(i+1)%8])
	                    transitions++;
	            }
	            if (transitions>=4) {                   //if neighbors contain more than one region, dilation ito this pixel is forbidden
	                table[item] = 0;
	            } else {
	            }
	        }
	        return table;
	    } // int[] makeFateTable	 
	 
	 /** dilate the UEP on one level by one pixel in the direction specified by step, i.e., set pixels to 255
	     * @param level the level of the EDM that should be processed (all other pixels are untouched)
	     * @param pass gives direction of dilation, see makeFateTable
	     * @param ip1 the EDM with the segmeted blobs successively getting set to 255
	     * @param ip2 a processor used as scratch storage, must hold the same data as ip1 upon entry.
	     *            This method ensures that ip2 equals ip1 upon return
	     * @param table the fateTable
	     * class variables used:
	     *  xCoordinate[], yCoordinate[]    list of pixel coordinates sorted by level (in sequence of y, x within each level)
	     *  levelStart[]                    offsets of the level in xCoordinate[], yCoordinate[]
	     * @return                          true if pixels have been changed
	     */
	private boolean processLevel(int level, int pass, Image8 ip1, Image8 ip2, int[] table,
	    int[] histogram, int[] levelStart, short[] xCoordinate, short[] yCoordinate) {
	        int rowSize = ip1.getWidth();
	        int height = ip1.getHeight();
	        int xmax = rowSize-1;
	        int ymax = ip1.getHeight()-1;
	        byte[] pixels1 = (byte[])ip1.getPixels();
	        byte[] pixels2 = (byte[])ip2.getPixels();
	        boolean changed = false;
	        for (int i=0; i<histogram[level]; i++) {
	            int coordOffset = levelStart[level] + i;
	            int x = xCoordinate[coordOffset];
	            int y = yCoordinate[coordOffset];
	            int offset = x + y*rowSize;
	            int index;
	            if ((pixels2[offset]&255)!=255) {
	                index = 0;
	                if (y>0 && (pixels2[offset-rowSize]&255)==255)
	                    index ^= 1;
	                if (x<xmax && y>0 && (pixels2[offset-rowSize+1]&255)==255)
	                    index ^= 2;
	                if (x<xmax && (pixels2[offset+1]&255)==255)
	                    index ^= 4;
	                if (x<xmax && y<ymax && (pixels2[offset+rowSize+1]&255)==255)
	                    index ^= 8;
	                if (y<ymax && (pixels2[offset+rowSize]&255)==255)
	                    index ^= 16;
	                if (x>0 && y<ymax && (pixels2[offset+rowSize-1]&255)==255)
	                    index ^= 32;
	                if (x>0 && (pixels2[offset-1]&255)==255)
	                    index ^= 64;
	                if (x>0 && y>0 && (pixels2[offset-rowSize-1]&255)==255)
	                    index ^= 128;
	                switch (pass) {
	                    case 0: if ((table[index]&1)==1) {
	                        pixels1[offset] = (byte)255;
	                        changed = true;
	                    }
	                    break;
	                    case 1: if ((table[index]&2)==2) {
	                        pixels1[offset] = (byte)255;
	                        changed = true;
	                    }
	                    break;
	                    case 2: if ((table[index]&4)==4) {
	                        pixels1[offset] = (byte)255;
	                        changed = true;
	                    }
	                    break;
	                    case 3: if ((table[index]&8)==8) {
	                        pixels1[offset] = (byte)255;
	                        changed = true;
	                    }
	                    break;
	                    case 4: if ((table[index]&16)==16) {
	                        pixels1[offset] = (byte)255;
	                        changed = true;
	                    }
	                    break;
	                    case 5: if ((table[index]&32)==32) {
	                        pixels1[offset] = (byte)255;
	                        changed = true;
	                    }
	                    break;
	                    case 6: if ((table[index]&64)==64) {
	                        pixels1[offset] = (byte)255;
	                        changed = true;
	                    }
	                    break;
	                    case 7: if ((table[index]&128)==128) {
	                        pixels1[offset] = (byte)255;
	                        changed = true;
	                    }
	                    break;
	                } // switch
	            } // if .. !=255
	        } // for pixel i
	        //IJ.write("level="+level+", pass="+pass+", changed="+changed);
	        if (changed)                    //make sure that ip2 is updated for the next time
	            System.arraycopy(pixels1, 0, pixels2, 0, rowSize*height);

	        return changed;
	    } //processLevel	 
	 
	 
	 /** Do watershed segmentation on a byte image, with the start points (maxima)
	     * set to 255 and the background set to 0. The image should not have any local maxima
	     * other than the marked ones. Local minima will lead to artifacts that can be removed
	     * later. On output, all particles will be set to 255, segmentation lines remain at their
	     * old value.
	     * @param ip  The byteProcessor containing the image, with size given by the class variables width and height
	     * @return    false if canceled by the user (note: can be cancelled only if called by "run" with a known ImagePlus)
	     */    
	private boolean watershedSegment(Image8 ip) {
	        int width = ip.getWidth();
	        int height = ip.getHeight();

	        byte[] pixels = (byte[])ip.getPixels();
	        // create arrays with the coordinates of all points between value 1 and 254
	        //This method, suggested by Stein Roervik (stein_at_kjemi-dot-unit-dot-no),
	        //greatly speeds up the watershed segmentation routine.
	        Histogram histogram = new Histogram(256);
	        ImageHelper.loadHistogram(ip, histogram);
	        int[] histogramValus = histogram.values;
	        int arraySize = width*height - histogramValus[0] - histogramValus[255];
	        short[] xCoordinate = new short[arraySize];
	        short[] yCoordinate = new short[arraySize];
	        int highestValue = 0;
	        int offset = 0;
	        int[] levelStart = new int[256];
	        for (int v=1; v<255; v++) {
	            levelStart[v] = offset;
	            offset += histogramValus[v];
	            if (histogramValus[v] > 0) highestValue = v;
	        }
	        int[] levelOffset = new int[highestValue + 1];
	        for (int y=0, i=0; y<height; y++) {
	            for (int x=0; x<width; x++, i++) {
	                int v = pixels[i]&255;
	                if (v>0 && v<255) {
	                    offset = levelStart[v] + levelOffset[v];
	                    xCoordinate[offset] = (short) x;
	                    yCoordinate[offset] = (short) y;
	                    levelOffset[v] ++;
	                }
	           } //for x
	        } //for y
	        // now do the segmentation, starting at the highest level and working down.
	        // At each level, dilate the particle, constrained to pixels whose values are
	        // at that level and also constrained to prevent features from merging.
	        int[] table = makeFateTable();

	        Image8 ip2 = ip.duplicate();
	        for (int level=highestValue; level>=1; level--) {
	            int idle = 1;
	            while(true) {                   // break the loop at any point after 8 idle processLevel calls
	                if (processLevel(level, 7, ip, ip2, table, histogramValus, levelStart, xCoordinate, yCoordinate))
	                    idle = 0;
	                if (idle++ >=8) break;
	                if (processLevel(level, 3, ip, ip2, table, histogramValus, levelStart, xCoordinate, yCoordinate))
	                    idle = 0;
	                if (idle++ >=8) break;
	                if (processLevel(level, 1, ip, ip2, table, histogramValus, levelStart, xCoordinate, yCoordinate))
	                    idle = 0;
	                if (idle++ >=8) break;
	                if (processLevel(level, 5, ip, ip2, table, histogramValus, levelStart, xCoordinate, yCoordinate))
	                    idle = 0;
	                if (idle++ >=8) break;
	                //IJ.write("diagonal only; level="+level+"; countW="+countW);
	                if (processLevel(level, 0, ip, ip2, table, histogramValus, levelStart, xCoordinate, yCoordinate))
	                    idle = 0;
	                if (idle++ >=8) break;
	                if (processLevel(level, 4, ip, ip2, table, histogramValus, levelStart, xCoordinate, yCoordinate))
	                    idle = 0;
	                if (idle++ >=8) break;
	                if (processLevel(level, 2, ip, ip2, table, histogramValus, levelStart, xCoordinate, yCoordinate))
	                    idle = 0;
	                if (idle++ >=8) break;
	                if (processLevel(level, 6, ip, ip2, table, histogramValus, levelStart, xCoordinate, yCoordinate))
	                    idle = 0;
	                if (idle++ >=8) break;
	                //IJ.write("All directions; level="+level+"; countW="+countW);
	            }
	        }
	        return true;
	    } // boolean watershedSegment
	 

	/** analyze the neighbors of a pixel (x, y) in a byte image; pixels <255 are
     * considered part of lines.
     * @param ip
     * @param x
     * @param y
     * @return  IS_LINE if pixel is part of a line, IS_DOT if a single dot
     */    
	int isLineOrDot(Image8 ip, int x, int y) {
        int result = 0;
        int width = ip.getWidth();
        int height = ip.getHeight();
        byte[] pixels = (byte[])ip.getPixels();
        int offset = x + y*width;
        int whiteNeighbors = 0;             //counts neighbors that are not part of a line
        int countTransitions = 0;
        boolean pixelSet;
        boolean prevPixelSet = true;
        boolean firstPixelSet = true;       //initialize to make the compiler happy
        for (int d=0; d<8; d++) {           //walk around the point and note every no-line->line transition
            if (isWithin(ip, x, y, d)) {
                pixelSet = (pixels[offset+dirOffset[d]]!=(byte)255);
                if (!pixelSet) whiteNeighbors++;
            } else {
                pixelSet = true;
            }
            if (pixelSet && !prevPixelSet)
                countTransitions ++;
            prevPixelSet = pixelSet;
            if (d==0)
                firstPixelSet = pixelSet;
        }
        if (firstPixelSet && !prevPixelSet)
            countTransitions ++;
            //if (x>=210&&x<=215 && y>=10 && y<=17)IJ.write("x,y="+x+","+y+": transitions="+countTransitions);
        if (countTransitions==1 && whiteNeighbors>=5)
            result = IS_LINE;
        else if (whiteNeighbors==8)
            result = IS_DOT;
        return result;
    } // int isLineOrDot	
	
	
	/** Delete single dots and lines ending somewhere within a segmented particle
     * Needed for post-processing watershed-segmented images that can have local minima
     * @param ip 8-bit image with background = 0, lines between 1 and 254 and segmented particles = 255
     */    
	void cleanupExtraLines(Image8 ip) {
        int width = ip.getWidth();
        int height = ip.getHeight();
        byte[] pixels =  (byte[])ip.getPixels();
        for (int y=0, i=0; y<height; y++) {
            for (int x=0; x<width; x++,i++) {
                int v = pixels[i]&255;
                if (v<255 && v>0) {
                    int type = isLineOrDot(ip, x, y);
                    if (type==IS_DOT) {
                        pixels[i] = (byte)255;                  //delete the point;
                    } else if (type==IS_LINE) {
                        int xEnd = x;
                        int yEnd = y;
                        boolean endFound = true;
                        while (endFound) {
                            pixels[xEnd + width*yEnd] = (byte)255;  //delete the point
                            //if(movie.getSize()<100)movie.addSlice("part-cleaned", ip.duplicate());
                            endFound = false;
                            for (int d=0; d<8; d++) {               //analyze neighbors of the point
                                if (isWithin(ip, xEnd, yEnd, d)) {
                                    v = pixels[xEnd + width*yEnd + dirOffset[d]]&255;
                                    //if(x>210&&x<215&&y==13)IJ.write("x,y start="+x+","+y+": look at="+xEnd+","+yEnd+"+dir "+d+": v="+v);
                                    if (v<255 && v>0 && isLineOrDot(ip, xEnd+dirXoffset[d], yEnd+dirYoffset[d])==IS_LINE) {
                                        xEnd += dirXoffset[d];
                                        yEnd += dirYoffset[d];
                                        endFound = true;
                                        break;
                                    }
                                }
                            } // for directions d
                        } // while (endFound)
                    } // if IS_LINE
                } // if v<255 && v>0
            } // for x
        } // for y
    } // void cleanupExtraLines	
	
    /** after watershed, set all pixels in the background and segmentation lines to 0
     */
	private void watershedPostProcess(Image8 ip) {
        //new ImagePlus("before postprocess",ip.duplicate()).show();
        byte[] pixels = (byte[])ip.getPixels();
        int size = ip.getWidth()*ip.getHeight();
        for (int i=0; i<size; i++) {
           if ((pixels[i]&255)<255)
                pixels[i] = (byte)0;
        }
        //new ImagePlus("after postprocess",ip.duplicate()).show();
    }	
	
	  /** Create an 8-bit image by scaling the pixel values of ip (background 0) and mark maximum areas as 255.
	    * For use as input for watershed segmentation
	    * @param ip         The original image that should be segmented
	    * @param typeP      Pixel types in ip
	    * @param isEDM      Whether ip is an Euclidian distance map
	    * @param globalMin  The minimum pixel value of ip
	    * @param globalMax  The maximum pixel value of ip
	    * @param threshold  Pixels of ip below this value (calibrated) are considered background. Ignored if ImageProcessor.NO_THRESHOLD
	    * @return           The 8-bit output image.
	    */   
	void make8bit(Image16 ip, Image8 typeP, boolean isEDM, float globalMin, float globalMax, double threshold, Image8 outMap) {
	        int width = ip.getWidth();
	        int height = ip.getHeight();
	        byte[] types = (byte[])typeP.getPixels();
	        double minValue = threshold;
	        double offset = minValue - (globalMax-minValue)*(1./253/2-1e-6); //everything above minValue should become >(byte)0
	        double factor = 253/(globalMax-minValue);
	        //IJ.write("min,max="+minValue+","+globalMax+"; offset, 1/factor="+offset+", "+(1./factor));
	        if (isEDM && factor >1./EDM.ONE) 
	            factor = 1./EDM.ONE;   // in EDM, no better resolution is needed
	        byte[] pixels = (byte[])outMap.getPixels();  //convert possibly calibrated image to byte without damaging threshold (setMinAndMax would kill threshold)
	        long v;
	        for (int y=0, i=0; y<height; y++) {
	            for (int x=0; x<width; x++, i++) {
	                v = Math.round((ip.getPixelValue(x, y)-offset)*factor);
	                if (v<0) pixels[i] = (byte)0;
	                else if (v<=255) pixels[i] = (byte)(v&255);
	                else pixels[i] = (byte) 255;
	            }
	        }

	        pixels = (byte[])outMap.getPixels();
	        //new ImagePlus("pre-threshold", outIp.duplicate()).show();
            for (int y=0; y<height; y++) {          //Set out-of-threshold to background (inactive)
                for (int x=0; x<width; x++) {
                    if (ip.getPixelValue(x,y) < threshold)
                        pixels[x+y*width] = (byte)0;
                }
            }
	        for (int i=0; i<width*height; i++) {
	            if ((pixels[i]&255) == 255)             //pixel value 255 is reserved
	                pixels[i]--;
	            if((types[i]&MAX_AREA)!=0)
	                pixels[i] = (byte)255;              //prepare watershed by setting "true" maxima+surroundings to 255
	        }
	    } // byteProcessor make8bit

	 public void findMaxima(Image16 inputMap, double tolerance, double threshold, Image8 maxMap, boolean isEDM)
	 {
	        //new ImagePlus("find maxima input", ip.duplicate()).show();
	        int width = inputMap.getWidth();
	        int height = inputMap.getHeight();

	        maxMap.clear();
	        
	        makeDirectionOffsets(inputMap.getWidth(), inputMap.getHeight());
	        globalMin = Float.MAX_VALUE;
	        globalMax = -Float.MAX_VALUE;
	        for (int y=0; y<height; y++) {         //find local minimum/maximum now
	            for (int x=0; x<width; x++) {      //ImageStatistics won't work if we have no ImagePlus
	                float v = inputMap.getPixelValue(x, y);
	                if (globalMin>v) globalMin = v;
	                if (globalMax<v) globalMax = v;
	            }
	        }
	        
	        maxPoints = getSortedMaxPoints(inputMap, maxMap, excludeEdgesNow, isEDM, globalMin, threshold); 
	        
	        analyzeAndMarkMaxima(inputMap, maxMap, maxPoints, excludeEdgesNow, isEDM, globalMin, tolerance);
	 }
	 
	 public void performWatershed(Image16 inputMap, double threshold, Image8 maxMap, boolean isEDM, Image8 outMap)
	 {
	        //new ImagePlus("Pixel types",typeP.duplicate()).show();
	        if (outputType==COUNT || outputType==POINT_SELECTION)
	            return;
	        
        
	        if (outputType == SEGMENTED) 
	        {                  
	        	//Segmentation required, convert to 8bit (also for 8-bit images, since the calibration may have a negative slope)
	            
	        	make8bit(inputMap, maxMap, isEDM, globalMin, globalMax, threshold, outMap);
	        	
	            //eliminate all the small maxima (i.e. those outside MAX_AREA)
	            cleanupMaxima(outMap, maxMap, maxPoints);    
	            
	          //do watershed segmentation
	            watershedSegment(outMap);              
	            
	            if (!isEDM) 
	            	cleanupExtraLines(outMap);       //eliminate lines due to local minima (none in EDM)
	            watershedPostProcess(outMap);                //levels to binary image

// MH: to do	            
//	            if (excludeOnEdges) 
//	            	deleteEdgeParticles(outMap, typeMap);
	            
	        } 
//	        else 
//	        {                                        //outputType other than SEGMENTED
//	            for (int i=0; i<width*height; i++)
//	                types[i] = (byte)(((types[i]&outputTypeMasks[outputType])!=0)?255:0);
//	            outIp = typeP;
//	        }
	        
//	        byte[] outPixels = (byte[])outMap.getPixels();
//	        //IJ.write("roi: "+roi.toString());
//	        if (roi!=null) {
//	            for (int y=0, i=0; y<outIp.getHeight(); y++) { //delete everything outside roi
//	                for (int x=0; x<outIp.getWidth(); x++, i++) {
//	                    if (x<roi.x || x>=roi.x+roi.getWidth() || y<roi.y || y>=roi.y+roi.getHeight()) outPixels[i] = (byte)0;
//	                    else if (mask !=null && (mask[x-roi.x + roi.getWidth()*(y-roi.y)]==0)) outPixels[i] = (byte)0;
//	                }
//	            }
//	        }
	    } 
	 
	 

}
